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ABSTRACT 

This  paper  constitutes  a  shorty  critical  survey  of 
the  current  applications  of  digital  computers  to  solving 
problems  in  structural  mechanics »   It  includes  a  brief 
description  of  typical  medium  and  large  size  computing 
systems  ( IBM-650  and  IBM=70i4.)j,  a  list  of  standard  mathema- 
tics problems  to  which  structural  mechanics  problems 
commonly  reduce  for  which  adequate  codes  exist  together 
with  a  list  of  these  codes  for  both  machines,  an  informal 
survey  of  mechanics  codes  completed  and  under  development 
at  representative  centers ^  a  detailed  discussion  of  three 
advanced  computer  applications  and  concludes  with 
recommendations  for  future  work  in  this  fielda 
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NYO-8670 
SOLVING  STRUCTURAL  MECHANICS  PROBLEMS 
ON  DIGITAL  COMPUTERS 

Introduction. 

The  purpose  of  this  paper  Is  to  give  a  critical  but 
abridged  survey  of  the  current  use  of  digital  computers 
for  the  solution  of  problems  In  structiiral  mechanics*   It 
was  not  possible  for  the  present  conference  to  produce  a 
complete  bibliographic  review  of  all  the  current  work  being 
done  In  this  area^  the  sheer  volume  Is  too  great  despite 
the  fact  that  the  first  electronic  automatic  computer,  the 
ENIAC  only  dates  back  to  19i|6«   At  the  same  time,  much  of 
this  work,  although  useful,  might  be  classified  as  routine. 
In  which  existing  mathematical  subroutines  are  strung 
together  to  effect  a  numerical  solution  of  a  problem  by  a 
method  which  is  a  more  or  less  direct  transcription  of 
standard  procedures  of  structural  mechanics  which  have  been 
used  in  the  past  as  the  basis  for  slide  rule  or  hand  computer 
calculations.   These  applications  will  merely  be  mentioned 
briefly  or  omitted  altogether.   The  only  applications 
considered  here  in  detail  are  of  a  more  advanced  character 
both  in  the  mechanics  and  in  the  mathematics. 

In  Section  A  there  is  a  brief  review  of  the  IBM-650  and 
TOi;  systems  as  typical  examples  of  intermediate  and  large  size 
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computers o  There  follows  a  list  of  mathematics  problems 
commonly  encountered  in  structural  mechanics  and  a  specific 
code  is  cited  for  each  of  these  problems  for  each  machine « 
The  purpose  of  this  is  to  furnish  a  comparison  of  capacities 
and  speeds,  to  indicate  the  present  state  of  the  art  in 
mathematical  subroutines  and  to  possibly  aid  the  reader  in 
selecting  a  code  for  a  problems 

The  remainder  of  Section  A  is  devoted  to  the  results  of 
an  Informal  survey  of  the  computer  applications  in  structural 
mechanics  being  carried  out  at  a  nuir>^'^r  of  universities, 
government  laboratories,  industrial  laboratories  and  engineer- 
ing consultant  finnSj,  with  emphasis  on  the  development  of 
general  codeso 

In  Section  B  three  particular  computer  applications  are 
studied  in  some  detail «   These  are 

(1)  Non-linear  deflections  of  shallow  spherical  shells,, 

(2)  Non-linear  bending  and  buckling  of  circular  plates, 

(3)  Plastic  torsion  of  cylinders» 

The  final  section,  contains  some  conclusions  and  suggestions 
for  future  work* 


A»   Machines.  Problems,  and  Programs. 

1»   Digital  computerse   It  is  appropriate  perhaps  to 
begin  with  a  brief  discussion  of  machines o   A  nxunber  of  makes 
and  models  are  in  current  use  j  however,  for  the  purposes  of 
this  paper  attention  will  be  focused  on  systems  of  two  sizes s 
the  large  size  as  typified  by  the  IBM-70lj.  and  the  inter- 
mediate size  as  typified  by  the  IBM-650.   The  nature  and 
capabilities  of  these  machines  is  assumed  to  be  sufficiently 
well-known  so  that  only  a  bare  recital  of  statistics  is 
neededo   Comparison  of  these  computers  on  specific  problems 
will  follow  later  and  will  sharpen  an  appreciation  of  their 
capacities  and  limitations. 

The  basic  IBM-b^O  possesses  a  magnetic  drum  memory  of 
2000  words o   Computations  are  performed  in  the  decimal  system 
and  there  are  10  digits  to  the  wordo   Average  access  time  to 
a  memory  location  is  2000  microseconds o 

The  IBM-70iJ^  has  a  magnetic  core  type  memory  usually  of 
8,192  words  but  also  is  available  with  I6g3^k   or  32  5768  words, 
each  36  binary  (approximately  10  decimal  digits)  bits  in 
lengtho   Addition  or  subtraction  requires  2l\.   microseconds 
and  multiplication  or  division  is  accomplished  in  2I4.O  micro- 
seconds. 

Programming  for  these  machines  can  be  carried  out  at 
either  of  two  levels.   At  the  level  of  the  more  expert  coder. 

For  a  complete  listing,  see  [1].   For  a  short  history 
of  digital  computers  see  [2]« 
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symbolic  code  language  Is  uaed,  which  itself  assembles  a 
machine  language  program  on  the  machine.   This  requires 
substantial  training  and  kn,owledge  of  the  machine o   Such 
codes  ares  for  example,  SAP  (SHA.RE  ASSEMBLY  PROGRAM)  for 
the  70l|.  and  SOAP  (  SYTXIBOLIG  OPTIMAL  ASSEM_BLY  PROGRAM)  for 
the  6S0b      At  the  level  of  the  inexperienced  machine  user 
who  wishes  to  devote  a  minimum  sjnount  of  time  in  learning 
to  code  „  there  are  automatic  ceding  systems  vjhich  more 
closely  resemble  the  language  of  mathematics  and  are 
designed  primarily  for  scientific  and  engineering  computa™ 
tion«   Among  these  are  IT  (INTERNAL  TRANSLATOR)  for  the 
650  and  FORTRAN  (FORMULA  TRANSLATOR)  for  the  70i|«   These 
codes  are  concise  to  write  and  themselves  generate  an 
efficient  symbolic  code  which  in  turn  generates  the  machine 
language  code.   There  is  a  so-called  FOR  TPJlNSIT  system 
which  makes  FORTRAN  available  to  650  users  by  translating 
FORTRAN  statements  Into  ITo 

2o  Mathematics  problems  and  programs ,  Problems  in 
structtiral  mechanics  frequently  reduce  to  one  or  another 
standard  mathematics  problem  of  a  type  for  which  adequate 
programs  already  exist  in  the  library  of  subroutines  for 
the  machine  in  use©  Although  the  author  wishes  In  large 
measure  to  discourage  the  "as  is"*  use  of  existing  codes, 
which  in  many  cases  would  be  inefficient  or  inadequate. 


particularly  for  larger  problems,  still  there  are  numerous 
problems  where  it  is  most  convenient  to  use  available 
programs  o 

Through  an  organization  of  users  such  as  SHARE  for  the 
Toil  and  the  IBM-650  Program  Library  codes  are  exchanged  and 
distributedo   The  libraries  of  codes  available  have  become 
so  voluminous  that  the  person  Interested  in  selecting  a 
program  for  solving  a  given  problem  needs  the  guidance 
furnished  by  the  experience  of  the  "old=hands"  at  his 
particular  computing  centero 

Typical  and  perhaps  most  common  among  the  standard 
mathematics  problems  to  be  solved  in  structural  mechanics 
for  which  good  programs  are  available  are  those  of 

(I)  Matrix  Inversion  and  solution  of  a  system 
of  linear  algebraic  equations 

(II)  Finding  the  eigenvalues  and  eigenvectors  of 
a  matrix 

(III)  Integrating  a  system  of  first  order  ordinary 
differential  equations 

(Iv)  Linear  programmingo 
Problems  Involving  the  solution  of  partial  differential 
equations  leading  through  difference  approximations  to  large 
systems  of  linear  (or  non-linear)  algebraic  equations  are 
common  today  but  in  general  one  must  select  a  method  and 
write  a  code  tailored  to  the  individual  problemo   Examples 
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of  these  will  be  given  In  Section  Bo 

In  the  following;,  specific  codes  are  cited  for  the 
650  and  70i|.  for  the  problems  listed  above o   This  is  done 
to  furnish  a  comparison  of  the  capacities  and  speeds  of 
the  machines,  to  indicate  something  of  the  present  state 
of  the  art.,  and  to  possibly  assist  the  reader  in  the 
selection  of  a  program  for  his  problem.   Although  these 
codes  are  considered  among  the  bestj,  there  are  others 
equally  good  both  for  these  arid  other  machines., 

(i)   Matrix  Inversion  and  simultaneous  linear 
equations  „   A  code  by  D..  V«  Sweeney  li|]  for  the  IBM-650 

inverts  matrices  of  order   <  1^.2     or  solves  b   sets  of 

p 
simultaneous  equations  with  n~  +  nb  <   17614.,   The  matrix 

elements  are  given  in  floating-point  fonrjo   Gaussian 

elimination  is  followed  and  coefficients  from  the  diagonal 

are  u''ed  regardless  of  size,  (which  can  cause  trouble  )o 

3 
The  matrix  inversion  is  carried  out  in  approximately  <.072  n 

seconds  plus  input -output  time.  Thus,  for  n  =  kOj  the 

time  required  for  the  calculations  is  over  1  l/k   hours o 

For  the  IBM-TOii,  an  ITY^J  program  [5.1  based  on  the 

Los  Alamos  LAS  665  code  available  through  SHARE  and  using 

Gaussian  elimination  will  calculate  the  inverse  of  a 

kO   X  kO   matrix  in  36  seconds  and  of  a  60  x  60  matrix  in 


For  a  general  sur-vey  of  methods  and  bibliography^  see 
[3], 
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2  minutes e   The  instructions  for  this  code  occupy  some 
225  words  of  memory  so  that  in  principle  the  order  of  the 
system  which  can  be  solved  is  limited  only  by  the  need  to 
store  coefficients  in  the  remainder  of  the  memory.   As  far 
as  storage  is  concerned,  therefore,  in  an  8,192  word  70l|., 
a  matrix  problem  of  order  89  could  be  attempted.   However, 
success  in  terms  of  accuracy  in  the  answer  is  dependent 
upon  the  condition  of  the  matrix.   Thus,  for  the  ill- 
conditioned  Hilbert  matrix  with  elements  h. .  =  (i+j-1)   , 
i,j  =  Ij  .«,,  n,   for  n  =  J4.  with  a  determinant  whose 
value  is   1,65  x  lO"'   about  5  significant  digits  can  be 
obtained  accurately  in  the  Inverse  matrix.   At   n  =  9 
with  a  determinant   9.72  x  10~^"^   in  value,  the  code  fails 
to  give  any  digits  accurately.   A  comparatively  well 
conditioned  matrix  such  as  \~~       with  elements 

r^.  =   i(n+l-j)/(n+l),   i  <  j,   Yji  =  -r^y      i.J  =  1,  •••»  n 

~2 
gives  for  n  =  60,   with  a  determinant  of   la61|  x  10 

about  5  significant  figures  accurately  in  the  time  of  two 
minutes  previously  quoted. 

(ii)  Eigenvalues  and  eigenvectors  of  a  real  symmetric 
matrix.   A  code  developed  at  Indiana  University  for  the  650 
(IBM-650  Library  Program  Abstract  not  yet  distributed)  based 
on  the  Jacobi  matrix  diagonalization  procedure  will  give  all 
eigenvalues  and  vectors  for   n  <  32   in  a  time  of  approxi- 
mately  (2o5  X  10"^ )n^  +  [|.  X  10"^  n^  minutes  and  eigenvalues 
-^  For  a  survey  of  methods  and  bibliography,  see  [6]  and  [3]» 
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only  for  32  <  n  <  Sh  ,    in  a  time  of  approximately   0<,006n-^ 
minutes  plu3  punch  out  time »   Thus „  to  calculate  all  elgen- 
values  and  vectors  for  n  =  30  would  take  about  3  hours  and 
20  minutes  « 

For  the  70l^.p  an  NYTJ  code  (NU  MXEV  -  available  through 
SHARE)  based  on  Given" s  method  [7!  will  handle  an  80  x  80 
matriXj,  given  8192  words  of  memory j,  and  a  120  x   120  mabrlx^, 
given  l6p38ij.  words  of  memory o   In  fixed  point  arithmetic 
(optional)  the  time  in  seconds  in  the  range   20  <  n  <  $0 
is  given  as   o3n  +  »075  n   to  give  both  eigenvalues  an4 
vectors o   For  n  =  30   the  time  would  be  approximately 
L  l/l\.   minutes  o   Loss  of  accuracy  with  increasing   n  is 
felt  first  in  the  vectors  then  in  the  eigenvalues  a  Far 
n  =  120j,   in  floating  pointp  double  precision  arithmetic 
would  probably  be  required  Just  to  get  the  largest  eigen- 
value o 

(iii)  Ordinary  differential  equations o   For  the  IBM»650, 
Ro  Wo  De  Sio  has  prepared  a  code  [8]  for  the  solution  of  a 
system  of  first  order  equations  j^   =  f  j  (Xpy-^P .,  00  py^^J  j, 
i  =  Ip  coop  Hj,   subject  to  given  initial  values  at  x  =  0, 
This  is  a  fixed-point  routine  based  on  v    predlctor-oorE'ector 
method  that  is  repeated  until  the  desired  convergence  is 
obtained  at  each  stepo   The  storage  required  is  380  locations 
plus  what  is  required  for  the  calculation  of  the  function  f^p 
which  in  turn  would  determine  the  running  time  per  x   incre- 
ment » 
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For  the  IBM»70i|.,  the  NYU  code  NU  DEQ  [9]  can  be  used 
to  integrate  the  above  systerrio   The  method  is  based  on  a 
variation  by  S„  Gill  ClO]  of  the  Runge-Kutta  procedur-eo 
This  is  a  floating-point  code  using  110  storage  locations 
for  the  instructions  exclusive  of  the  calculation  of  the  f., 

(Iv)  Linear  programming o  This  is  a  newer  and  less 
standard  type  of  problem  than  the  preceding j,  but  one  which 
is  rapidly  growing  in  importance o  Ws  shall  refer  to  a 
problem  of  this  type  later  as  well  as  a  problem  in  quadratic 
programming  occurring  in  structural  mechanics 9  The  linear 
programming  problem  in  general  is  that  of  minimizing  or 
maximizing  a  linear  expression  in  a  set  of  unknowns 
X,  jooopX  J,   subject  to  a  number  m  of  side  conditions  which 
consist  of  either  linear  algebraic  equations  or  inequalities 
involving  the  unknowns o  Unrestricted  variables  (as  to  sign) 
can  always  be  replaced  by  non-negative  variables  by  intro- 
ducing new  unknowns o   Inequalities  similarly  can  always  be 
replaced  by  equations  by  increasing  the  number  of  unknowns,, 
The  codes  cited  next  assume  these  changes  to  be  madeo 

For  the  IBM^650,  there  is  a  code  by  L„  So  Woo  [11] 
based  on  the  Modified  Simplex  Method  El2]  which  handles  a 
maximi-im  of  97  equations  in  an  unlimited  number  of  variables  o 
The  time  required  is  stated  to  vary  from  l\.   to  13  minutes 
per  iteration  during  the  first  [|.0  Iterations » 

For  the  IBM-YOi;  a  code  by  Orchard- Hayes  „  Judd  and 
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Cutler,  available  through  SHARE,  also  based  on  the  Modified 
Simplex  Method,  will  handle  up  to  2^5  equations  in  any 
niiinber  of  variables »   A  general  estimate  of  time  required 
is  not  available!  however,  it  is  stated  that  an  example 
involving  ^0   equations  and  100  variables  required  20 
minutes  for  solution, 

3«  Structural  mechanics  problems.   In  an  effort  to 
ascertain  what  use  is  currently  being  made  of  digital 
computers  for  the  solution  of  structural  mechanics  problems, 
the  author  conducted  an  informal  survey  of  a  number  of 
university,  government  and  industrial  sites  where  computing 
is  done  on  these  problems.   The  list  included  Brown  University, 
M,I»To,  University  of  Illinois,  University  of  Houston,  New 
York  University,  National  Bureau  of  Standards,  David  Taylor 
Model  Basin,  Aberdeen  Proving  Grounds,  Boeing  Airplane 
Company,  Gulf  Oil  Corporation,  IBM  Corporation,  General 
Electric  Research  Laboratories,  General  Electric  Plight 
Propulsion  Laboratory  Department,  Westinghouse  Research 
Laboratories,  M.  W,  Kellogg  Company,  American  Bridge 
Division  of  U.  S.  Steel,  Gibbs  and  Cox,  and  the  firm  of 
Howard,  Needles,  Tammen  and  Bergendoffo   This  list  could  have 
been  extended  indefinitely  but  is  considered  to  be  a  fair 
sample. 

Before  reviewing  the  technical  information  obtained 
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in  the  course  of  this  survey „  it  may  be  appropriate  to 
consider  the  need  for  exchange  of  information  of  this  kindo 
The  number  of  centers  and  the  number  of  individuals  Involved 
in  writing  programs  for  the  express  purpose  of  solving 
problems  in  structural  mechanics  is  multiplying  so  rapidly 
that  any  static  attempt  to  compile  a  bibliography  or  list- 
ing of  codes  would  be  of  value  for  only  a  short  time  and 
would  likely  be  outdated  before  it  appeared  In  print o 
Accordingly  it  seems  clear  that  what  is  needed  is  an 
organization  of  interested  groups  to  provide  a  mechanism 
for  the  continuing  interchange  of  codes  and  other  informa- 
tion pertaining  to  the  solution  of  structural  mechanics 
problems  on  digital  computers o   This  could  be  accomplished 
through  the  publication  by  such  an  organization  of  a  periodic 
newsletter 5,  for  example^  containing  reports  sent  in  on  a 
voluntary  basis  by  members  on  their  codes  and  activities© 
Representatives  of  the  cooperating  groups  could  meet 
periodically  in  conjunction  with  the  regular  meetings  of 
an  appropriate  technical  society  such  as  the  Applied 
Mechanics  Division  of  the  AoSoMoEo   Precisely  this  sort  of 
organization  has  been  set  upj,  for  example,  for  the  exchange 
of  nuclear  reactor  codes  and  is  known  as  the  Nuclear  Codes 
Group  (NCG)o   Its  newsletter  is  distributed  through  one  of 
the  member  groups j,  the  AEG  Computing  and  Applied  Mathematics 
Center  at  New  York  University.   This  organization  of  users 
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has  been  assisted  by  IBM  which  acts  to  distribute  decks  of 
cards  as  requested  by  members  for  actual  codes  for  IBM 
machines  which  have  been  submitted  to  the  NCG* 

Returning  now  to  the  survey  of  structural  mechanics 
codeSj,  it  is  perhaps  most  instructive  for  the  purposes  of 
the  present  discussion  to  describe  the  situation  in  general 
terms,  as  it  appeared^,  rather  than  to  attempt  to  list 
specific  codes o 

Obviously  every  paper  that  appears  in  a  journal  giving 
numerical  results  for  a  specific  problem  and  acknowledging 
the  use  of  a  computer  in  itself  constitutes  a  computer 
application,,   However j,  one  may  make  a  general  distinction 
between  one-shot  applications j,  or  codes  that  serve  only 
oncej,  and  general  codes  that  are  designed  to  cover  a  whole 
class  of  problems.   The  first  type  is  noteworthy  when  a 
new  method  is  developed  or  the  application  is  unusualo 
It  is  the  second  type  of  codej,  however,  that  is  of  the 
most  interest  to  us  here.   It  appears  that  people  are  now 
in  the  first  phases  of  writing  such  general  piorpose  codes 
for  structural  mechanics* 

Starting  with  the  mathematically  simplest  problems, 
we  find  that  there  are  codes  for  the  elastic  stress  analysis 
of  two  and  three  dimensional  trusses  and  rigid  frames^  under 

^  Extensive  work  on  computer  analysis  of  such  structures 
has  been  carried  out  by  Lives I'^y  in  England »  See,  for 
example,  [13]  and  also  bibliography  in  [173e 
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static  loads  (University  of  Houston,  Gulf  Oil  Corporation, 
American  Bridge  Division  of  U.  S.  Steel  among  others). 
Related  to  these  are  codes  for  the  analysis  of  thermal 
expansion  stresses  in  arbitrary  three  dimensional  piping 
systems  (Ma  W.  Kellogg  Company  among  others).   There  is  a 
code  for  the  elastic-plastic  analysis  of  trusses  calculat- 
ing both  stresses  and  displacements  up  to  the  plastic  limit 
load  (New  York  University),   Ordinary  plastic  limit  analysis 
of  trusses  and  framed  structures  can  be  carried  through  using 
standard  linear  programming  codes  [li;],   A  code  has  been 
developed  for  the  minimiim  weight  design  of  frames  (Brown 
University  [l5])o   Codes  for  the  elastic-plastic  response 
to  air  blast  loads  of  such  structures  as  beams,  circular 
plates,  rectangular  membranes  have  been  prepared  at  the 
Ballistic  Research  Laboratory  at  Aberdeen  Proving  O-roundSo 

A  number  of  codes  for  the  solution  of  structural  problems 
of  especial  concern  to  the  civil  engineer  have  been  developed 
at  the  University  of  Illinois  for  the  ILLIAC  under  the  general 
direction  of  N.  M,  Newmark,   As  reported  on  by  Veletsos  [16], 
these  include  the  stress  analysis  of  simply  supported  skew 
I-beam  bridges,  determination  of  natural  frequencies  and  modes 
for  the  vibration  of  multi-story  building  frames  with  flexible 
girders,  dynamic  response  of  building  frames  under  wind 
pressures,  blast  loads  and  ground  disturbances,  and  the 
dynamic  response  of  simple  span  highway  bridges  under  the 
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passage  of  heavy  vehicles^,  among  others.   Bridge  analysis 
programs  have  also  been  developed  at  the  University  of 
Houstono   In  general  the  civil  engineering  field  has  been 
active  in  the  development  and  use  of  computer  codes  for 
analysis  and  design  for  several  years  »   As  an  indication 
of  this  activity  we  cite  the  May  1956  issue  of  Civil 
Engineering  which  is  entirely  devoted  to  digital  computers 
and  their  applications,.   According  to  this  source,  over  100 
civil  engineering  firms  and  l\.0    state  highway  departments 
have  their  own  computers o   A  Committee  on  Electronic  Computa- 
tion of  the  i^tructural  Division  of  the  AoSoGoEo  is  currently 
conducting  a  survey  of  computer  codes » 

In  the  aircraft  industry  there  are  general  codes  for 
flutter  analysis  and  the  determination  of  the  natiiral 
frequencies  and  modes  of  delta  wings  and  "built-up"  structures 
formed  from  plates,  bars  and  stringers  (Boeing  Aircraft  among 
others ). 

The  analysis  of  cylindrical  shells  and  shells  of 
revolution  is  of  widespread  interest  in  such  connections 
as  the  design  of  submarines,  missiles,  turbines,  reactors, 
A  code  for  the  calculation  of  the  strength  of  stiffened 
cylindrical  shells  was  developed  at  the  David  Taylor  Model 

Basin.   At  the  General  Electric  Gas  Turbine  Division  several 

r-~- ■ — ^ — ~ 

Ro   Wo    Clough   [17]    has  made    a  survey  of   structural   analysis 
applications  of  computers   with  a  good  bibliography. 
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codes  have  been  developed  for  the  small  deflection  theory 
of  shells  under  axially  symmetric  and  also  ijnsymmetrlc 
loadings.   Of  considerable  interest  is  a  G.Eo  code  which 
will  reportedly  analyze  a  built  up  structure  of  three  shells 

of  revolution  joined  in  series,,   The  heart  of  the  GeEo  shell 

th 
programs  is  a  routine  to  solve  systems  of  8   order  linear 

ordinary  differential  equations© 

Methods  and  codes  for  the  solution  of  non-linear  bending 
and  buckling  problems  in  elastic  plates  and  shells  have  been 
developed  at  New  York  University  and  M.IoT.   Some  of  these 
developments  will  be  reviewed  later  in  this  paper  as  special 
examples,, 

Codes  have  been  and  are  being  developed  for  the  solution 
of  1|   order  biharmonic  type  boundary  value  problems  (NYU  and 
M.Iol'o)  and  problems  associated  with  2   order  quasi-linear 
partial  differential  equations  (NYU)o   The  former  have 
application  in  two-dimensional  elasticity  problems  including 
plane  strain,  plane  stress  and  plate  problemso   The  latter 
has  various  applications  in  mechanics,  for  example  in  plastic 
torsion  with  strain-hardening,  which  will  be  reviewed  later 
as  a  special  example » 
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B.   Some  Case  Histories o 

The  birds-eye  view  of  the  survey  of  the  last  section 
is  of  necessity  lacking  in  details   To  compensate  for  this, 
the  present  section  is  a  report  on  three  specific  computer 
applications  with  which  the  author  is  familiare   All  are 
currently  being  carried  out  at  the  AoEaCo  Computing  and 
Applied  Mathematics  Center  at  New  York  University,   It  is 
believed  that  these  applications  are  interesting,  first, 
because  they  are  concerned  with  problems  of  a  difficult  and 
advanced  character,  both  mechanically  and  mathematically 
and,  second,  because  of  the  light  they  shed  on  the  various 
roles  which  the  computer  can  play  in  research, 

1»   Non-linear  deflections  of  shallow  spherical  shells  <> 
The  mechanical  problem  is  the  determination  of  the  deflections 
of  a  thin  spherical  cap  under  uniform  external  pressure  and 
clamped  around  the  edge,,   To  account  for  experimentally 
observed  loads  below  those  predicted  by  the  linear  theory  of 
buckling  and  the  snap-through  form  of  buckling  mode   the 
problem  must  be  treated  as  one  in  the  non-linear  theory  of 
buckling  following  von  Karman ,  Tsien  and  FriedrichSo   In  [l8] 
the  problem  is  reduced  to  solving  the  following  pair  of 

non-linear  second  order  ordinary  differential  equations 

T — — — — — ■  ■  — — 

For  a  complete  report  see  [18]  including  references 
cited  in  the  present  discussion* 
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L  a  =  p  {ay  +  Pe^) 

(1) 

L  Y  =  p  (©^  -  a^) 

where   ^=®de9de®   ^^  ^^^   differential  operator,   ©, 
the  ratio  of  tx.o  polar  angle  to  the  total  angular  opening 
is  the  independent  variable,   P  is  a  loading  parameter, 
p   is  a  geometric  parameter  proportional  to  the  ratio  of 
the  shell  rise  to  thickness,   a  =  a(6)   is  proportional  to 
the  slope  of  the  shell,   y  -  y(®)   is  a  stress  function  and 
all  quantities  are  dimensionlesse   Deflections  and  stresses 
are  obtained  from  a   and   y»   The  boundary  conditions  are 
a,  =  Y  =  0   at  the  center,   ©  =  0,   and  a  =  1,   ^  -  v  y  -  0 
at  the  edge,   ©  =  1,   where   v   is  Poisson's  ratio. 

Although  solutions  for  this  problem  had  been  obtained 
in  the  past  by  perturbation  and  by  power  series  methods, 
they  unfortunately  covered  only  a  small  range  of  the  parameters. 
The  problem,  therefore,  was  to  devise  a  method  of  numerical 
solution  for  machine  calculation  which  would  cover  a  wide 
range  of  parameters,  from  very  shallow  shells  (near  plates) 
to  shells  with  a  high  rise,  from  low  loads  to  buckling  loads 
(in  some  way  to  be  determined)  to  perhaps  post-buckling  loads 
and  the  possible  numerical  determination  of  non-adjacent 
equilibrium  configurations.   In  other  words  what  was  desired 
was  a  code  which  would  enable  one  to  study  the  behavior  of  a 
whole  range  of  shells  under  widely  varying  loads, 
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A  power  series  raethod  of  solution  \jas  attempted  first  and 
carried  through  as  follows.   Inserting  expansions  of  the  form 

00-1  OO        0-,  1 

a  =  X^  ^r  ®     *       '''   ^ ^n  ®  ^"^ 

11=1   '  n=l 

into  (1),  recurrence  relaticns  were  obtained  for  the   a 

and   Y   in  terms  of  the  initial  coefficients   a^   and   y,  • 
'n  1       '1 

The  boiondary  conditions  at  the  edge  reduced  to 

I__  a   -  1  =  0  ,         1_  f2n  -  (H-v)y^]  =  0    (3) 
n=l  ^  n=l  ^ 

which  were  in  effect  then  simultaneous  equations  for   a 

and   Y-i   fo^  given  values  of  the  parameters   p  and   P<, 

For  fixed  values  of   p  and   P,   the  numerical  procedure 

consisted  in  at  first  choosing  trial  values  for  a^   and 

Yn  9   calculating  successive  coefficients   a   and   Yv^   i^ 

the  series  from  the  recurrence  relations  until  a  convergence 

criterion  was  satisfied,  and  ther.  evaluating  the  functions  in 

(3)  to  see  whether  these  equations  were  satisfied  to  a 

required  degree  of  accuracy.   Improved  values  of  a,   and 

Y,   were  obtained  systematically  by  a  Newton-Raphson  type 

of  iteration  toward  the  solution  of  the  equations  (3)» 

li^Then  these  equations  were  satisfied  the  displacements  and 

stresses  were  computed  from  the  power  series  expressions 

evaluated  from  the  final  values  of  a   and  y   obtained. 

n        n 

___________ 

Solutions  by  power  series  on  a  computer  have  also  been 
obtained  for  this  problem  at  M.I.T,  by  Welnitschke  [193 
under  Reissner. 
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For  a  given  shell   (p  value),  the  computations  were 
begun  with  a  small  load   (P  value)  for  which  the  deformed 
state  was  nearly  spherical  so  that  a  good  initial  estimate 
of  the   cl^jY    values  could  be  made.   The  Newton-Raphson 
process  then  rapidly  gave  the  roots   a..   and  y^      of  (3)  to 
the  desired  accuracy.   For  successively  higher  values  of   P^ 
the  initial  guess  for  a..   and  y,   was  found  by  extrapola- 
tion of  roots  obtained  for  preceding   P  values.   This  was 
aided  by  a  plot  of  points  in  an  a..  ,y-. -plane  for  increasing 
P  values  for  each  shello 

The  calculations  were  carried  out  on  the  AEG  UNIVAC  at 
New  York  University  using  a  code  of  some  300  words  in  length. 
It  did  not  prove  practical  to  have  the  initial   a^ ,y   guess- 
ing programmed  due  to  extreme  sensitivity  to  correct   a,  jYt 
choices  especially  for  larger  values  of   p   and   Po   Human 
judgement  was  needed  and  close  cooperation  between  man  and 
machine  was  required  for  good  results. 

For  a  given  shell,  the  load  was  increased  in  small 

Increments  and  solutions  obtained  according  to  the  described 

procedure.   The  loading  process  terminated  by  a  lack  of 

convergence  in  one  of  two  ways.   Either  the  Newton-Raphson 

process  failed  to  converge  (convergence  to  an  accuracy  of 

1  X  10~^  in  (3)  was  usually  obtained  in  under  10  iterations 

if  at  all)  or  else  the  higher  coefficients   a  .v   calculated 

&  n  'n 

from  ass-umed   o-i  sYi   values  failed  to  converge  (i.e«  to 
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become  less  than  1  x  10°°^  in  100  terms  )o 

Typically  the  process  ended  suddenlyo   Good  convergence 
was  obtained  for  a  succession  of  Increinents  In   P  for  a 
given  shell o   Then  for  any  slight  increase  in   P  no  solution 
could  be  obtained  by  the  process  through  a  selection  of 
0-,  sy-i      values  in  the  neighborhood  of  one  extrapolated  from 
values  for  slightly  lower  values  of  P<.   It  was  as  evident 
as  possible  in  a  numerical  process  that  no  neighboring 
solution  existed  to  the  preceding  oneSo   The  buckling  load 
was  then  assumed  to  lie  between  the  highest  attained  value 
of  P  for  which  the  process  converged  and  the  next  larger 
value  attemptedo   This  was  completely  analogous  to  the 
situation  in  a  laboratory  test  when  the  buckling  load  is 
reached  aad  a  new  equilibrium  configuration „  substantially 
different  from  the  previous  ones,,  is  attained^   As  examples 
of  this  numerical  "buckling"  behaviors  for   p  ~   32   and 
P  =  lo25  the  procedure  converged  to  a  solution  in  3  iterations* 
Increasing  P  to  lo3»  no  convergence  in  the  Newton-Raphson 
process  was  in  evidence  after  10  iterationso   The  buckling 

load   P    for  p  =  32  was  therefore  assumed  to  lie  between 

cr       r   ^ 

P  =  lo25  and  P  =  lo3<.  For   p  =  26ol  and  P  =  lo55  conver- 
gence occxirred  in  2  iterations  and  only  ^0  terms  of  the  power 
series  were  needed „  while  for  P  =  lo6     neighboring  ^.^.^''^l 
values  could  not  be  found  to  lead  to  convergent  power  series 
even  after  99  terms »   The  buckling  load  in  this  case  was 
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assiuned  to  lie  between  P  -  l»b>5  and  P  ~   I060  The  results 
were  in  good  agreement  with  experiments  done  by  Kaplan  and 
Fling  [20 ]o   As  a  time  estimate  for  solving  a  given  shell 
under  a  given  load^,  the  time  per  Iteration  was  roughly  one 
minute  in  the  case  where  50  terms  in  the  power  series  were 
neededo   Thus  in  the  cases  cited,  2  or  3  minutes  of  computa- 
tion led  to  a  complete  listing  on  tape  of  the  radial 
deflection  and  all  stress  components  for  the  top  and  bottom 
surface  of  the  shell  as  functions  of  angle  for  the  given 
loade 

Interestingly  enough,  it  was  possible  for  some  values 
of  p  to  pick  up  solutions  corresponding  to  loads  in  excess 
of  the  buckling  load.   This  was  not  done  by  continuously 
following  the  solutions  as  the  loads  are  Increased  from 
below  the  buckling  loado   As  explained  above  this  was  not 
posslbleo   However,  one  could  locate  the  post-buckling 
branch  of  the  curve  in  the  ci-ipYt   plane  by  taking  advantage 
of  the  two  parameter  character  of  the  problem,,   For  low 
enough  values  of  pj,   the  shell  approximates  a  flat  plate 
and  no  buckling  occurs  so  that  one  may  continuously  locate 
solutions  up  to  any  load,  including  those  in  excess  of 
buckling  loads  for  shells  of  larger  p  value,,   Having 
attained  a  solution  for  small   p   and  large   P,   one  can 
then  vary   p  keeping   P  constant  and  continuously  proceed 
to  the  solution  for  values  of   p  for  which  the  attained  P 
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is  in  the  post-buckling  range©   One  can  then  decrease   P 
for  fixed  p  and  find  a  sequence  of  solutions  for  the 
loads  between  the  buckling  load  and  attained  loado  Thus, 
numerically  one  could  locate  different  branches  of  solutions, 
say  in  the   o,,  jY-i   plane,  and  from  these,  different  branches 
of  typical  load- deformation  cvirves  (say  for  the  center  of 
the  spherical  cap)  corresponding  to  different  ranges  of 
loadso 

The  computations  clearly  showed  certain  effects  which 
had  not  been  predicted  theoretically  before  and  for  which, 
although  experimentally  observed^  there  existed  little 
precise  data.   These  effects  were  the  dependence  of  buckling 
mode  on   p^   the  shape  parameter  together  with  an  associated 
non-monotone  dependence  of  buckling  load  on   p   and  the 
existence  of  a  boundary  layer  for  higher  po   Thus,  it  was 
foiindj,  from  plots  of  the  radial  deflection,  that  for  p  <  23 
approximately,  the  ordinary  mode  of  buckling  occurred  in 
which  the  maximiim  deflection  is  at  the  center  and  falls  to 
zero  at  the  edgeo   Increasing   p   still  more  leads  to 
flattening  of  this  maximxam  and  transition  to  a  second  mode 
where  the  raaximxam  moves  from  the  center  progressively  toward 
the  edge,  and  where  rapid  changes  in  deflection  and  stress 
define  a  boxindary  layer o   This  mode  seems  to  persist  until 
around   p  =  52  where  a  third  mode  begins  to  take  shape  with 
two  maxima,  one  at  the  center.   Plots  of  buckling  load  against 
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p  revealed  a  peaking  behavior  with  the  peaks  occurring 

at  the  values  of   p  where  the  transitions  in  mode  occurred. 

Generally  speaking,  the  power  series  solution  was  no 
longer  satisfactory  when  it  came  to  large  values  of   p   or 
large   P   (post-buckling),  due  to  the  necessity  for  increas- 
ingly large  numbers  of  terms  in  the  expansions.   It  was  also 
dependent  upon  finding  close  initial  estimates  for   a^   and 
Y-|   and  could  proceed  therefore  only  by  continuous  variation 
of  parameters  and  solutlonso   The  next  example  discusses  a 
method  of  "iterative  differences"  which  was  applied  subse- 
quently to  the  present  shell  problem  and  proved  to  be  much 
faster  and  relatively  insensitive  to  initial  guessese   Using 
this  technique  solutions  for  a  larger  range  of   p  and   P 
were  obtained, 

2»   Non-linear  bending  and  buckling  of  circular  plates 
[21,22  3o   A  thin  circular  elastic  plate  is  considered, 
subjected  to  either  a  uniform  lateral  pressure  or  a  uniform 
radial  edge  thrust  with  clamped  or  simply  supported  edges » 
The  von  Karman  plate  equations  in  polar  coordinates  for 
rotatlonally  symmetric  deformations,  reduce  to 


2 
a(x)  y(x)    "    P  X 


ik) 


J    -1    J 

where   L£x-r---r-x   is  the  second  order  differential 
dx  X  dx 
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operator  as  previously,   x   is  the  dlmensionless  radius, 
a(x)   is  a  dimenslonless  slope,   Y^^)   is  the  dimenslonless 
radial  derivative  of  the  Airy  Stress  Function,  and  P  is 
the  dimenslonless  lateral  load.   Various  boundary  conditions 
are  considered,  in  particular  we  may  here  cite  the  bending 
problem  (A)  with  clamped  edge,,  for  which 

ail)   =0  ,    [f|»  vrl^^,  -0  , 

and  the  buckling  problem  (E)  with  clamped  edge,  for  which 

a(l)  =  0  ,  y(1)  =  P  >  0  o 

In  both  cases  by  sjrmmetry  and  regularity,   a  =  y  ~  0  at 
the  center,   x  =  0« 

The  numerical  method  of  solution  consisted  in  approxi- 
mating the  equations  (I4.)  by  finite  differences,,   The  resulting 
system  of  non-linear  algebraic  equations  for  the  unknown 
function  at  the  radial  mesh  points  x.  =  ih,   i  =  1,  •<>•,  m-1 , 
were  then  solved  by  the  following  "interpolated"  iteration 
process.  Bo\indary  conditions  are  Imposed  at  the  end  points 

X   =  0 ,   X  =  1 « 
o     •    m 

Lett...g  L  stand  for  the  difference  operator  correspond- 
ing to   L  a  simple  iteration  process  is  defined  by 

L  a^(x^)  =  -  a^_^  (X.)  y^,^  (x^)  -  P  x^ 

^   P  ,0<i<m   (^) 

LYjx,)=|af  (X,) 
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starting  from  an  Initial  guess  for  a  (x. ),   y  (x. )   the 

Ox        O    1 

first  of  equations  (5)  defines  a  system  of  in-2  linear 

equations  for  m-2  unknowns   a  (x.  ) ,  Prom  the  a-(x.)si 

the   Y-1  (x.  )   can  be  found  from  the  second  of  equations  (5) 

as  the  solution  of  a  linear  system.   However,  analysis  and 

calculations  show  that  this  process  will  converge  only  for  a 

small  range  of  the  parameters   P   or   P  and  will  definitely 

diverge  for  values  of  P  and   P  above  certain  predicted 

values  o 

An  improved  procedure  (which  can  be  motivated  on 

mechanical  grounds)  is  obtained  if  one  interpolates  at  the 

n   step  between  the  previous  iterate  a  , (x. )   and  the  one 

n-l   1 

found  as  the  solution  of  the  first  set  of  equations  in  (5), 
and  then  uses  the  interpolated  function  as  the  new  iterate. 
Thus,  one  sets 

a^(x^)  =  cja^Cxj.)  +  n  -UJ)   ^^.^(x^)        (6) 

where   a   denotes  the  solution  of  the  equations 
n 

L  a^(x.)  =  -  a^^^(x^)  y^^:^U^)    -  P  xf  . 

Here  (^    is  the  interpolation  parameter,  which  lies  between 
0  and   1 o 

The  best  value  to  use  for  u)    is  determined  empirically 
by  setting  up  a  convergence  criterion  such  as 


max 
0<i<m 
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axid  finding  the  number  of  iterations  N(cj)   required  as  a 

function  of   a)  ,   Typically  one  finds  a  graph  with  a  sharp 

and  sensitive  minimuin  which  is  a  function  of  the  parameter 

value  of   P   or  P.   This  defines  the  optimum  Ci3  value e 

Fortunately  the  optimum  iJ    was  found  to  be  relatively 

insensitive  to  the  mesh  size  so  that  using  a  coarse  mesh 

and  little  computing  time  one  could  obtain  by  machine 

experiments  the  proper  LJ    values  to  be  usedo 

In  problem  (A),  the  initial  iterate  was  taken  to  be 

the  converged  solution  for  the  next  smaller  load  value. 

In  problem  (E)  there  is  the  possibility  of  more  than  one 

solution  and  the  method  should  be  applied  with  this  in 

mindo   For   F  <  F    (where   F    is  the  lowest  eigenvalue 

o  o 

of   the    linear  buckling   problem)    the    only   solution   is   that 

of  direct  compression,      a(x.)    =  0,      yix^)   =  Px^o      The  value 

of     F   ,      however,      is  not  needed   a  priori*      One   can  start 

with  a  non-zero  estimate   for     a(x    )g      for      F  <  P        the 

iteration  will   converge   to  the   trivial  solution     a(x.)    =  0, 

By  increasing     P      in  small   increments   from  zero,   at  some 

P^     the   iteration  will  converge   to  a  non-trivial  solution 
o  ° 

and     F        may  be   taken  as   an   estimate   of     F    .      From  here    on 
as     P     is    increased  above     P   ,      the   previous   non-trivial 
solution  is   used  as    the   first   estimateo     When    the    iterations 
converge   for  a  given  value    of    the   load,    the   stresses  and 
deflections   are  computed, 
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The  calculations  were  carried  out  on  the  AEC  UNIVAG 
at  Nevi  York  University,  using  10,  20  and  $0   mesh  points. 
For  50  points  some  500  words  of  memory  were  required  for 
the  complete  code.   The  method  was  found  to  give  excellent 
convergence  for  a  practically  unlimited  range  of  loads »   In 
case  (A),  values  of  the  parameter   P  as  high  as   P  =  7OOO 
were  easily  attainedo 

For  comparison  the  method  of  power  series  previously 
discussed  was  also  used  for  the  present  problem  in  case  (A)., 
The  solutions  obtained  agreed  closely  with  those  of  the 
present  method.   It  was  however  found  that  the  power  series 
method  failed  for  much  lower  loads,  around   P  =  2000 «   More- 
over, the  difference-iteration  method  was  faster  in  general 
by  a  factor  of  ten  (thus  most  solutions  were  obtained  in 
under  a  minute)  and  was  more  automatic  being  relatively 
insensitive  to  the  initial  guesses,  which  could  therefore 
be  programmed. 

The  results  gave  a  clear  definition  of  the  boundary 
layer  effect  first  discussed  by  Friedrichs   and  Stoker 
and  the  completeness  of  the  data  made  possible  an  empirical 
determination  of  the  width  of  the  boundary  layer  as  a 
function  of  the  loading  parameter* 

3,   3»   Plastic  torsion  of  cylinders »   The  problem  considered 
is  that  of  torsion  of  a  prismatic  cylinder  of  a  material  with 
a  well-defined  yield  such  as  mild  steelo   The  work  to  be 
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summarized  here,  due  to  \i o  S.  Dorn  and  the  author p  is  still 
in  progress  and  will  be  presented  elsewhere  with  a  complete 
discussion  of  method  and  results. 

We  give  a  brief  review  of  the  formulation  of  the 

o 

problem   o      Taking   the      x,y-plane    in   the   cross-section  and 
the      z-axis      parallel   to  the  generators    of   the   cylinder   the 
St.   Venant   assumptions  for    the   displacement   components 
u,   V,    w      are 

u  =  -yzO   ,      V   =  xzQ   ,        w   =  w(x,y,e) 

where   ©  is  the  angle  of  twist  per  unit  length  of  the 
cylindero   The  strains  are  zero  except  for  the  shear 
strains  given  by 

The  components  f  ,  TT  of  the  shear  stress  can  be  obtained 
from  a  stress  function  t(^»7ii®^   ^5 

^.  =  if  '  -y  =  -  il  '«' 

and   equilibrium  is    thereby  automatically  satisfied   for  any 
function     \j(     vanishing   on   the   boundary   of    the    cross-section 
(assuming  a   simply  connected  region). 

Making   the   assumption   of   an   elastic-perf ectly  plastic 

p : 

For  a  more  detailed  discussion  of  the  plastic  torsion 
problem  see  Prager  [23]© 
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material,  Hooke's  law  provides  the  stress-strain  relation 
as  long  as  the  condition   Igrad  \|r|  =  J (^~)      +  ^ji^     ^  ^      ^^ 
satisfied  locally,  where   k   is  the  yield  stress  in  pure 
shear.   Yielding  takes  place  at  points  where  the  condition 
Igrad  \|(  I  =  k   is  satisfied.   As   6   increases  from  zero  the 
solution  is  everywhere  elastic  until  the  yield  condition  is 
satisfied  locally,  following  which  plastic  zones  grow  tmtil 
ultimately  (for  an  oo  angle  of  twist)  the  entire  cross- 
section  would  become  plastiCo 

The  mathematical  problem  is  to  determine  for  each  © 
a  function   \lf(x,y,9)   vanishing  on  the  boundary,  continuous 
with  continuous  first  partial  derivatives  such  that  every- 
where  Igrad  \1(  I  <  k   and  such  that  wherever   jgrad  \1;|  <  k, 
the  elastic  compatibility  equation  d^^/d^^^'   +  ^^^/ij^   =  ' -2G^© 
is  satisfied  where   G    is  the  elastic  modulus  of  rigidity. 

The  problem  of  solution  is  made  difficult  by  the  -unknown 
boundary  between  the  elastic  and  plastic  regions  and  the 
distinct  partial  differential  equations  to  be  satisfied  in 
these  regions.   No  satisfactory  systematic  niomerical  procedure 
has  been  devised  to  handle  this  problem.   The  difficulties 
arise  because  of  the  assumption  of  perfect  plasticityo  This 
assximption,  which  represents  an  approximation  to  the 
phenomenon  of  yield  observed  in  metals  like  mild  steel,  is 
introduced  to  simplify  tne  mathematical  theory  in  those 
problems  where  the  entire  material  is  yielding  or  the  extent 
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of  the  plastic  zone  Is  known  or  satirsfles  certain  a  priori 
conditions  such  as  those  of  limit  analyslSa   Howeverp  In 
cases  like  the  present  one  the  assiunptlon  of  perfect 
plasticity  leads  to  unnecessary  mathematical  complications. 

To  eliminate  these  difficulties j,  one  may  use  a  plasti- 
city theory  for  a  material  >diose  stress-strain  curve  Is  of 
the  continuous  transition  variety <>   These  curves  have  an 
initial  elastic  slope  and  bend  continuously  approaching 
either  the  yield  stress   k  of  perfect  plasticity  as  ,an 
asymptote  or  increase  beyond  this  to  incorporate  the  effects 
of  strain  hardeningo 

Ramberg  and  Osgood  [2l4,J  introduced  a  set  of  stress- 
strain  curves  depending  on  three  parameters  which  fit  a 
variety  of  metals  and  include  work-hardening  effect So 
For  simple  shear,  the  stress  f  and  strain  y   (in  the 
present  notation)   are  related  by  an  equation  of  the  form 

Y=.l.  [1+^]^  (9) 

where  n  Is  a  positive  integer  parameter,   G   and  k  are 
essentially  the  elastic  rigidity  and  yield  stress©   For 

fixed  n,  kj,   G  ^   the  stress  is  a  monotone  increasing 

o 

function  of  strain  which  flattens  as  it  reaches  k  in  value 
but  Increases  beyond  k   in  value  as  the  strain  Increase So 
As   n  — >  00   the  stress-strain  curve  approaches  the  broken 
line,  flat  yield  cixrve  of  a  naterial  which  behaves  elastlcally 
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with  shear  modulus   G   up  to  the  shear  stress   k  and  then 

o 

yields   under  constant   shear   stress      k»      The   larger     n     is, 
the  more   nearly   is  this   relation  approxiraatedo 

To   appropriately  generalize    (9)   to    the  case  of  several 
stress-components   as    in  torsions    one  may  either  proceed  to 
formulate  a  deformation-type    (finite)   stress-strain  law 
following  Nadai    [25]    or  a  flow-type    (differential)    stress- 
strain  law  following   Prager   [23,    25,    26]    and  H,    J.    Laning 

9 

(unpublished  paper)© 

Appropriate  relations  of  the  deformatioi'i  type  are 
given  by 

(t^+r|)^ 

^x   2a^  ^=^     j^2n   ^   X 

(^  2jn  (10) 

^.  =  ^'^*^^=^.     • 

In  terms  of  the  reduced  stress  function  ^   =  \|;/k,   the 
use  of  (10)  leads  to  the  following  partial  differential 
equation  for  <^    I 

(l+l2n^2ni?/2(n-l)j^^^^^^^^j^^g2(n-l)^^^^  ^ 
+  (l+l2n+2nU?yl2(n-l)j^^^  ^  ^  ^ 

■    ■   ■■   '   rii  t   ■   .   I   ■  t  .■.iii.Jii.a     I,.  I  ■   ■   III  .i—i.i-.    II .   ■   ■   ■   1   I   I   I       ■   L       ■  I         1  !■»   1   »■■  ■  ^w 

The  Ramberg-Osgood  curves  were  also  used  by  Huth  [2?] 
in  a  comparison  of  flow  and  deformation  theories  in 
plastic  torsion*   His  calculations  were  carried  out 
for  a  square  cross-section  using  a  low  n  value  and 
a  relatively  coarse  mesho 
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(11) 


— 2        — 2  ""2 

where      S      =   tl/j      +  ^t        and   a  coinina  followed  by      x     or     y 

indicates   partial   differentiation  with  respect   to    that 

variable.      For   a   simply  connected  cross-section   this 

equation   is    to   be    solved   subject   to    the    condition     -^  =  0 

on   the  boundary.      The  numerical    solution  of   this   problem 

will  be  considered  after   the    formulation  which  follows  of 

the   torsion  problem  in  terms   of   a   flow  theory. 

Appropriate    stress-strain  relations  of    the   flow  type 

are   given  by 

P      P   n-1 

2G^Y    =  ^    +  (2n+i)      ^  J   —  (r  r    +  r  r  )  t 

o^         y  jj2n  '"-x  •-  X  y     y'  ''y 

where   the   dot   denotes    differentiation  with  respect    to      9 
which   is    assumed  to    increase  monotonically  and   so   plays 
the    role   of  the   time   variable.      Using   these  relations,    the 
partial   differential    equation  for      ■!ir(x,y,0)     becomes 

[l+^?^s2(^"lJ(2n+l)]i^^   +   [2(n+im2(n-l);j^^j^^j   ^^^^  ^ 
+   [l+1??yS2('^~l)(2n+l)]i,^   +   rs2(^-l>(l?,^f,^^+^,^^,^y)    + 
+  Ki,^](2n+1)?,^  +   [S^^''""^^(^,^"iF,^y+i;',yi,yy)   +  K^ , ^ ]  ( 2 H +1 ) ^ , ^ 


G 

__o_ 

k 


,  (13) 


where      K  =  2(n-l  )3^^''"^  ^[lFf^)|f,^^+2iV,3^t,yi,^+^?yiF,yy]  o 
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This    equation   is   subject    to    the  boundary  conditions     ^(xs,y,@)-0 

for   a  simply  connected  cross-section  and   for     6  >  0» 

The   deformation  theory  equation  for     "^     is   of  the  form 

2G  0 
^'xx  ^'xy  ^*yy  k 

where   A,  B,  C   are  functions  of  "^,   and  xjr,   and  is  a 

X       y 

quasi-linear  second  order  partial  aifferential  equation  of 
elliptic  type.   A  numerical  procedure  was  devised  at  NYLT 
by  S»  Schechter  for  handling  general  equations  of  this  type 
and  coded  by  R.  Wernick  for  the  IBM-TOi;,   This  code  was 
successfully  used  for  the  present  problem.   The  following 
is  a  brief  description  of  the  procedure,  the  code  and  the 
results  obtained. 

Substituting  differences  for  the  derivatives  in  the 
equation  one  obtains  a  set  of  non-linear  algebraic  equations 
for  the  values  of  ^     defined  on  the  points  of  a  rectangular 
mesh.   As  in  the  previous  example,  an  iterative  procedure  was 
used  for  the  solution  of  these  equations.   Starting  from  an 
initial  guess   \|(^°'   at  each  point,  the  new  value  at  a  given 
mesh  point  was  calculated  as  follows*   The  coefficients  A, 
B,  C   of  the  difference  equations  are  evaluated  from  the 
given  ^^°   values.   The  differences  for  the  second  order 
terms  are  calculated  from  the   \1;     values  as  well  except 
for  the  value  at  the  given  mesh  point  which  is  calculated 
so  as  to  satisfy  the  equation.   This  value   \[r   constitutes 
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the  uncorrected  new  Iterate  value  at  the  point o  To  obtain 

=-(1) 
the  new  iterate  value   t,     at  ^'h©  points,  again  a  weighted 

average  of  the  form 


is  taken  where  the  parameter  /UJ  riow  is  foimd  to  give  the 
best  convergence  for  value?  greater  than  one,  furnishing 
a  so-called  extrapolated  value.   Moving  to  the  next  mesh 
point  proceeding  in  a  prescribed  order.,  say  row  by  roWp 

the  coefficients   A,  Be,  C   and  second  order  differences 

—  (o ) 
are  next  computed  from  the  old  values   (|f     together  with 

"(1) 
the  new  value  of  \jr     calculated  at  the  preceding  point* 

In  this  way  the  latest  Iterate  values  are  always  used  and 

each  step  furnishes  an  explicit  value  for  the  next  Iterate 

at  the  given  mesh  pointo 

=—(  n+l  1  ^~f  n  1 

For  a  convergence  criterion ^   f       and  ■^  are 

compared  at  each  point  of  the  mesh  and  the  maximum  of 

1'^^     "^    I   taken  over  all  the  points  is  taken  as  a 

st 
measure  of  the  error  at  the   (n+l)    iteration  or  "sweep" 

of  the  mesho   In  the  present  process  the  calculation  was 
stopped  when  this  quantity  was  less  than  1  x  lO"  » 

To  gauge  the  size  of  the  computational  problem  corres- 
ponding to  a  two  dimensional  problem  of  this  complexity s, 
the  following  data  is  informative «   The  Schechter-Wernlck 
code  was  designed  for  more  general  equations  in  which  the 
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coefficients  A,  By  C  could  be  functions  of  ^ j,    ts^s  ts  » 
X  and  y  and  also  for  boundary  conditions  involving  the 
normal  derivative  of  ^      as  well  as   v|fo   It  also  allowed 
for  multiply-connected  domains  and  permitted  variable  mesh 
spacing  from  point  to  point.   The  basic  code  for  all  this 
took  up  some  2800  words  of  memory.   In  addition,  some  l500 
words  were  required  for  general  purpose  subroutines  such  as 
inputs  output  and  edite   Over  l+OO  words  were  used  in  our 
problem  in  routines  for  the  calculation  of  the  specific 
coefficient  functions  kg   Bj,  C,,   and  for  a   31+  x  31+  point 
mesh  (the  limit  of  the  code)  some  3i4-00  words  of  storage 
were  required.   Thus,  virtually  the  entire  8,196  word 
memory  of  a  standard  70i+  was  required. 

The  method  was.  checkedout  against  the  case  of  a  circular'  cross- 
section  where  the  solution  could  be  found  analytically  and 
was  found  to  give  excellent  agreement »   The  calculations  were 
then  carried  out  for  the  case  of  a  square  cross-section  where 
no  analytic  results  are  available  but  the  character  of  the 
solution  for  the  limiting  case   (n  — >  oo  )   of  an  elastic- 
perfectly  plastic  material  is  well  known  from  the  soap-film, 
sand-hill  analogy  of  Nadai,   According  to  this  analogy,  as 
0,   the  angle  of  twist,  is  increasedj,  the  plastic  zones  will 
develop  in  from  the  centers  of  the  sides  of  the  squares  and 
gradually  fill  out  the  four  triangles  formed  by  the  diagonals 
of  the  square.   The  material  in  strips  adjacent  to  the 
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diagonals  constitutes  a  residual  elastic  region.   In  the 
plastic  region,  the  stress  intensity  measured  by  S   should 
remain  constant  equal  to   ko 

The  calculations  bore  this  out  and  indeed  for   n  =  1? 
with  2G  9/k  =16,   among  the  highest  values  of  these 
quantities  used,  the  plastic  zones  were  sharply  defined 
regions  of  essentially  constant   So 

With  a  20  X  20  mesh  and  n  =  13^   starting  from  an 
arbitrary  initial  guess  of  ^    =  l/8  at  all  points, 
convergence  to   1  x  10~    required  about  50  or  60  sweeps 
at  i    seconds  per  sweep  and  so  cook  about  6  minutes  for  a  given 
©  value.  With  n  =  1?   (so  that  the  coefficients  A,  B,  C 
involved  3h  powers  of  the  differences  for  ■ij;,    and   \j(,  , 

which  were  calculated  by  multiplications  rather  than  logar- 
ithms) the  time  went  up  to  about  10  seconds  per  sweep  for  a 
total  time  of  about  10  minutes  for  a  whole  calculation« 
Starting  from  a  previous  solution  for  a  smaller  ©  value 
for  the  initial  guess  ^    ,   the  time  was  cut  by  about 
one -half  in  each  case. 

For  large  values  of  2G  ©/k,   above  16,  the  elastic 
zone  along  the  diagonals  became  too  small  for  adequate 
definition  using  a   20  x  20  mesh  in  the  square  and 
a  314.  X  3^  mesh  was  used,  which  was  at  the  limit  of  the 
given  code.   With  this  mesh,  values  of   20  ©/k   as  high  as 
25  were  reached.   However,  with  this  nximber  of  points  the 
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time  per  sweep  went  up,  for  n  =  17j)   to  about  25  seconds 
and  solution  time  for  a  given  angle  of  twist  to  20  to  30 
minutes,  revealing  the  enormous  amount  of  calculation 
involved. 

It  Is  Intended  to  compare  the  results  of  this  solution 
of  the  deformation  theory  equation  with  a  solution  of  the 
flow  equation  (13)  and  to  further  compare  predictions  of 
the  two  types  of  theories  for  other  cross-sections,.   Such 
comparisons,  followed  by  further  comparisons  in  two- 
dimensional  problems,  it  is  hoped  will  lead  to  the  selection 
of  a  valid  stress-strain  relation  covering  through  parameters 
a  number  of  metals,  which  can  be  used  as  a  basis  for  complete 
elastic-plastic  analyses  of  continuous  structures  using 
computers. 

The  flow  equation  (13)  considered  as  one  In  the  three 
independent  variables  x,  y   and  G   is  of  third  order  and 
is  essentially  parabolic   Equations  of  this  type  In  which 
the  "time"  derivative  is  not  given  explicitly  but  appears 
implicitly  as  the  solution  of  a  space  variable  elliptic  type 
problem  is  encountered  also  In  meteorology  problems.   The 
numerical  problem  of  solution  is  a  formidable  one  and  there 
is  as  yet  comparatively  little  known  theoretically  or 
empirically  about  methods  of  solution.   A  method  which 
suggests  Itself  Is  to  begin  v/lth  a  small   ©  value,  use 
the  elastic  solution  to  calculate  the  coefficients  of  the 
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highest  order  derivatives  and  then  solve  the  linear  2nd 
order  problem  for  \l;  ,   These  values  can  be  used  in  a  time 
step  to  advance  the  value  of  ^     to  the  next  0  value  and 
then  repeat  the  process.   Probably  accuracy  will  be  improved 
by  some  form  of  averaging  of  time  step  values.   It  does  not 
appear  that  there  is  any  numerical  stability  problem  here  as 
there  is  in  time  dependent  problems  of  the  heat-flow  type. 
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C.   Conclusions  and  Beginnings, 

Any  number  of  self-evident  things  can  be  said  about 
the  future  importance  of  computers  in  mechanics  or  any 
other  science.   To  avoid  saying  these  again  is  difficult. 
However,  the  facts  speak  for  themselves.   The  number  of 
computers  available  is  already  large  and  growing  constantly. 
The  number  of  people  to  whom  programming  and  coding  and 
operating  a  computer  is  routine  is  growing  at  a  pace  where 
in  a  short  time  it  will  be  as  familiar,  if  not  quite  as 
portable,  sji  aid  to  the  student  as  the  slide  rule.   The 
speed  and  capacity  of  computers  is  being  enlarged  (on  the 
order  of  a  factor  of  10)  with  each  generation  of  computers 
(on  the  order  of  5  years)  to  the  point  where  it  is  safe  to 
assume  that  with  a  practical  amount  of  analysis,  time  and 
money  we  can  soon  solve  any  of  the  mechanics  problems  we 
currently  can  formulate.   Automatic  programming  routines 
will  allow  the  user  to  solve  his  problems  even  on  the 
largest  computers  with  relative  ease  and  a  minimum  of  effort. 

What  is  needed  is  intelligent  use  of  these  machines. 
As  has  been  indicated  there  is  a  place  for  problem  solving 
by  simple  application  of  existing  mathematical  subroutines 
appropriately  combined.  There  is  also  a  place,  although  a 
disappearing  one,  for  simple  adaptation  of  methods  and  formulas 
of  an  older  era  in  structural  mechanics  when  approximation 
was  dictated  by  the  computational  limitations  of  human  mind 
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or  haxido   The  challenge  to  be  met  by  today's  expert  Is  to 
forinulate  his  problem  and  his  method  of  solution  to  the 
limits  of  accuracy  of  existing  computers  and  ultimately  to 
formulate  as  complete  and  acciirate  a  description  of  the 
problem  as  Is  possible  based  on  existing  physical  knowledge,. 

By  virtue  of  the  accuracy  and  completeness  of  the 
results  obtainable  to  problems  properly  formulated  and 
solved  on  computers  revolutionary  new  procedures  become 
feasibleo   For  example,  the  computer  may  be  used  to  experi- 
ment on  a  structure  in  varieties  of  ways  and  with  a  control 
of  variables  never  possible  in  a  physical  laboratoryo   If 
the  model  programmed  into  the  computer  has  been  conceived 
in  broad  enough  terms  to  include  the  essential  parameters 
so  as  to  fully  describe  the  processes  involved^  every 
conceivable  combination  of  loads  j,  material  properties^  etCoj 
can  be  investigated  and  with  a  speed  and  simplicity  other- 
wise impossible.   The  example  of  the  spherical  shell  caps 
cited  previously  provided  just  such  an  opportunity o   Shells 
were  picked  and  loaded  to  buckling  and  beyond  as  they 
suggested  themselves  to  the  experimenter  almost  as  would 
have  been  done  in  an  actual  laboratory,  giving  measurements 
of  stresses  and  deflection  impossible  to  obtain  on  actual 
specimens. 

Again,  the  computer  may  be  used  to  select  a  design  by 
comparing  the  behavior  of  a  variety  of  structures  of  a  many 
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parameter  famllyo  An  outstanding  application  of  this 
approach  today  lies  in  the  field  of  design  of  turbine 
blades, 

A  basic  application  of  the  computer  will  be,  as  in 
the  plastic  torsion  problem^  to  compare  the  predictions 
of  competing  theories  in  critical  experiments. 

In  all  of  this  it  is  clear  that  new  methods  of 
analysis  must  be  developed!   mathematical  procedures  must 
be  suited  to  the  mechanical  problems  concerned  and  designed 
for  efficient  machine  usee 

At  the  present  time,  the  most  standard  tool  of 
numerical  analysis  in  problems  in  continuxom  mechanics  is 
the  difference  equation  approximation  to  the  differential 
equation.   This  leads  to  a  system  of  simultaneous  algebraic 
equations,  whichj,  even  in  two  dimensions  for  a  single 
function  for  say  a   10  x  10  mesh  is  of  order  100   and  in 
three  dimensions  for  a   10  x  10  x  10  mesh  is  of  order  lOOOo 
It  is  true  that  even  in  the  non-linear  case  5  iterative 
procedures  of  solution  (such  as  described  for  the  plastic 
torsion  problem)  exist  wherein  the  computational  complexity 
at  each  step  (say  in  going  from  point  to  point)  is  independent 
of  the  total  nvimber  of  unknowns »   However,  the  stability  of 
the  procedure  and  the  accuracy  of  the  results  becomes  increas= 
Ingly  doubtful  as  the  order  growso   One  is,  therefore ^  faced 
with  the  dilemma  that  by  putting  more  points  into  the  mesh  to 
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Improve  in  principle  the  accuracy  of  one's  solutlonj,  one 
loses  accuracy  for  reasons  of  the  accumulation  of  errors 
In  the  number  of  numerical  operations  required  for  solution 
using  the  given  method  of  solution. 

Ideally  what  is  needed  are  methods  of  solution  in  which 
greater  accxiracy  is  obtained  as  one  Increases  the  amount  of 
arithmetic  in  calculating  the  solution.  This  is  true  of  the 
Monte  Carlo  methods  [28]  which  simulate  the  physical  processes 
and  wherein  the  accuracy  of  the  results  increases  as  more 
trials  are  run  or  more  games  are  playedo   For  processes  like 
diffusion  where  a  physical  model  can  be  set  up  describing  the 
probabilistic  behavior  of  Individual  particles  the  Monte 
Carlo  method  is  most  suitable  and  is  being  used  successfullyo 
Unfortunately,,  the  deterministic  (or  what  we  today  consider 
deterministic)  problems  of  solid  mechanics  do  not  in  their 
physical  models  suggest  natural  Monte  Carlo  procedureso   One 
might  indeed  devise  Monte  Carlo  procedures  for  the  solution 
of  the  specific  differential  equations  of  a  given  mechanics 
problemj  but  this  is  not  appealing  as  it  runs  counter  to  the 
spirit  of  the  method  which  is  used  partly  to  avoid  the 
differential  equation  description.   Still  It  may  be  that  a 
more  basic  microscopic  formulation  of  mechanics  problems  in 
terms  of  fundamental  physical  laws  will  one  day  be  possible 
and  that  at  that  time  the  Monte  Carlo  methods  will  perhaps 
be  the  natural  ones. 
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Lacking  Monte  Carlo  procedures  and  limited  in  scope 
by  difference  equations  there  are  few  guideposts  to  methods 
for  solving  complex  structures  in  three  dimensionso   One 
suggestion  which  the  author  should  like  to  make^  and  whose 
practicality  has  yet  to  be  tested ^  is  a  return  to  the  direct 
methods  of  the  calculus  of  variations  ■=-  in  modern  dress© 
Traditionally,  for  linear  problems,  this  has  meant  approxi- 
mating the  solution  by  a  linear  combination  of  m  selected 
functions  and  then  minimizing  the  appropriate  energy 
expression  with  respect  to  the  coefficients  of  these 
functions  leading  to  a  system  of  linear  equations  of  order 
m«   The  result  is  the  best  combination  of  the  selected 
functions. 

To  describe  an  alternative  procedure  let  us  consider 
the  specific  example  of  finding  the  displacements   Uj  v,  w 
throughout  an  elastic  solid  under  prescribed  surface 
tractions  and  displacements*   We  replace  the  continuous 
body  by  a  rectangular  lattice  of  n   points  and  seek  the 
displacements  at  these  points ,  a  total  of   3n  unknowns 
(this  would  be  an  enormous  nxrmberj   3000  in  the  case  of 
even  a  modest   10  x  10  x  10  mesh  of  points )»   Next,  the 
potential  energy  is  expressed  in  terms  of  these  unknowns,, 
Using  first  differences  for  strains  and  combinations  of 
these  for  stresses  this  yields  a  quadratic  expression  in 
3n  unknowns  to  be  minlmizedo   The  boundary  conditions 
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can  be  used  to  reduce  the  number  of  unknowns  in  this 
expression  or  can  be  left  (particularly  when  the  conditions 
are  on  stresses)  as  subsidiary  linear  conditions  on   some 
of  the  unknowns o 

If  one  next  proceeded  to  use  the  calculus  to  find  the 
conditions  for  a  minimum,  one  would  be  led  to  a  system  of 
3n  linear  equations j  more  or  less,  for  the  unknowns.   These 
would,  of  course,  comprise  a  system  of  approximating  differ- 
ence equations  for  the  partial  differential  equations 
governing  the  displacement  functions   u,  Vj,  w<,   In  fact, 
this  is  recognized  as  a  useful  and  consistent  way  of 
deriving  such  equations.   However^  this  procedure  leads  to 
the  problem  of  solving  a  linear  system  of  large  order  which 
is  what  we  wish  to  avoid. 

It  is  suggested  that  instead  of  the  above  procedure, 
attention  be  concentrated  directly  on  the  problem  of 
minimizing  the  quadratic  function  in  the  unknowns.   Any 
linear  elasticity  problem  can  of  course  always  be  reduced 
to  this  form.   The  problem  of  finding  direct  methods  of 
minimizing  a  quadratic  function  subject  to  linear  constraints 
is  one  that  is  currently  receiving  considerable  attention 
under  the  name  of  quadratic  programming  [29-35] o   Various 
algorithms  have  been  devised  for  proceeding  step-wise  to  a 
solution  of  the  problem.   The  most  recent  and  promising  of 
these  due  to  Wolfe  [353  is  based  on  a  simple  modification  of 
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the  simplex  method  of  linear  programmingo  These  methods^ 
however  J,  are  intended  for  problems  wherein  the  variables 
are  subject  to  linear  inequalities  and  hence,  which  caraiot 
be  reduced  via  the  calculus  to  the  solution  of  linear 
equations o   The  programming  techniques  are  usually  deemed 
unnecessary  where  no  inequalities  are  present  as  is  the 
case  in  the  elasticity  problemo   However,  the  algorithm 
devised  by  Wolfe  can  be  applied  to  this  case  in  principle 
to  find  after  a  finite  number  of  iterations  the  set  of 
values  for  Uj,  Vp  w  which  minimizes  the  potential  energy o 
How  this  procedure  compares  in  time  and  accuracy  with  the 
solution  of  linear  equations  procedure  remains  to  be 
studiedo  At  the  present  time  the  IBM<=70l|  code  prepared 
by  Wolfe  for  the  quadratic  programming  problem  and  based 
on  the  previously  cited  code  of  Orchard-Hayes p  Judd  and 
Cutler,  limits  the  sum  of  unknowns  and  subsidiary  condi- 
tions to  255o   Thusj,  without  modification  it  would  only 
be  feasible  to  study  at  most  two  dimensional  elasticity 
problems  and  even  here  the  limitations  are  stringent  since 
the  procedure  would  more  than  double  the  nixmber  of  physical 
unknowns  to  be  determined,   ( However g  for  structures  like 
frames  the  number  of  unknowns  allowed  would  seem  ampleo) 

Clearly  this  is  a  case  where  the  existing  mathematical 
procedure  should  be  studied  in  the  light  of  the  given 
mechanics  problem*   There  is  no  question  that  larger  machines 
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will  raise  the  boxinds  on  the  niunber  of  linknowns  in  the 
method  as  it  stands  to  practical  levels »   On  the  other 
hand  J,  by  tailoring  the  method  to  the  problem  there  is  a 
good  chance  that  the  matrices  can  be  reduced  in  size  to 
the  point  where  good  use  of  present  machines  can  be 
made  with  essentially  this  type  of  solution  of  elasticity 
problemso 
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